from osgeo import gdal_array

# 使用gdal 库加载图片到numpy库
src = gdal_array.LoadFile("thermal.tif")  # 把遥感影像作为一个数组读取
# 根据类别数目将直方图分割成20个颜色区域
classes = gdal_array.numpy.histogram(src, bins=20)[1]
# 颜色查找表的记录必须为len(classes) +1
# 声明RGB元组
lut = [[255, 0, 0], [191, 48, 48], [166, 0, 0], [255, 64, 64], [255, 115, 115], [255, 116, 0], [191, 113, 48], [255, 178, 115], [0, 153, 153], [29, 115, 115], [0, 99, 99], [166, 75, 0], [0, 204, 0], [51, 204, 204], [255, 150, 64], [92, 204, 204], [38, 153, 38], [0, 133, 0], [57, 230, 57], [103, 230, 103], [184, 138, 0]]
# 1111111111111111111111111111111
# 分类初始值
start = 1
# 创建一个RGB颜色的JPEG 输出图片（3个维度 行列 ）
rgb = gdal_array.numpy.zeros((3, src.shape[0], src.shape[1], ), gdal_array.numpy.float32)

# 处理所有类并声明颜色
for i in range(len(classes)):
# logical_and (numpy 用法 逻辑与或非)
    mask = gdal_array.numpy.logical_and(start <= src, src <= classes[i])
    for j in range(len(lut[i])):  # len(lut[i]) == 3
        rgb[j] = gdal_array.numpy.choose(mask, (rgb[j], lut[i][j]))
    start = classes[i] +1
output = gdal_array.SaveArray(rgb.astype(gdal_array.numpy.uint8), "classfied.jpg", format="JPEG")